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Perturbative schemes utilizing a spectral moment expansion are well known and extensively used 
for investigating the physics of model Hamiltonians and real material systems (in combination with 
density functional theory). However, such methods are not always reliable in various parameter 
regimes such as in the proximity of phase transitions or for strong couplings. Nevertheless, the 
advantages they offer, in terms of being computationally inexpensive, with real frequency output at 
zero and finite temperatures, compensate for their deficiencies and offer a quick, qualitative analysis 
of the system behavior. In this work, we have developed such a method, that can be classified 
as a multi-orbital iterative perturbation theory (MO-IPT) to study N-fold degenerate and non 
degenerate Anderson impurity models. As applications of the solver, we have combined the method 
with dynamical mean field theory to explore lattice models like the single orbital Hubbard model, 
covalent band insulator and the multi-orbital Hubbard model for density-density type interactions in 
different parameter regimes. The Hund’s coupling effects in case of multiple orbitals is also studied. 

The limitations and quality of results are gauged through extensive comparison with data from the 
numerically exact continuous time quantum Monte Carlo method (CTQMC). In the case of the 
single orbital Hubbard model, covalent band insulators and non degenerate multi-orbital Hubbard 
models, we obtained an excellent agreement between the Matsubara self-energies of MO-IPT and 
CTQMC. But for the degenerate multi-orbital Hubbard model, we observe that the agreement with 
CTQMC results gets better as we move away from particle-hole symmetry. We have integrated 
MO-IPT with density functional theory based electronic structure methods to study real material 
systems. As a test case, we have studied the classic, strongly correlated electronic material, SrVOa. 

A comparison of density of states and photo emission spectrum (PES) with results obtained from 
different impurity solvers and experiments yields good agreement. 

PACS numbers: 71.27.+a, 71.10.Fd, 71.10.-w, 71.30.-|-h, 71.15.Mb 


I. Introduction 

The development of efficient methods to solve quan¬ 
tum impurity problems, especially those involving multi¬ 
ple orbitals, has been a significant research direction in 
the field of theoretical condensed matter physics. Sub¬ 
sequent to the development of the dynamical mean field 
theory (DMFTjP, which is exact in the limit of infinite 
dimensions and an excellent local approximation in fi¬ 
nite dimensions, the importance of obtaining reliable so¬ 
lutions to general quantum impurity problems has in¬ 
creased further. 

Within the DMFT framework, a lattice model may 
be mapped onto a quantum impurity embedded in a 
self-consistently determined host. The impurity problem 
may then be solved by a variety of techniques including- 
numerically exact methods like quantum Monte Carlo 
(QMC), numerical renormalization group (NRG), exact 
diagonalization (ED) and density matrix renormaliza¬ 
tion group (DMRG) or semi-analytical methods like iter¬ 
ated perturbation theory (IPT), local moment approach 


(LMA), non-crossing approximation (NCA) and fluctua¬ 
tion exchange approximation (FLEX). Each method has 
its own advantages as well as pitfalls. For example, 
QMCP is a numerically exact method, but is computa¬ 
tionally expensive. It yields data on the Matsubara axis 
(or imaginary time) so to obtain dynamical quantities 
such as the density of states and transport quantities, 
analytic continuation of the data to real frequencies is 
essentiaP, which is a mathematically ill-posed problem. 
Additionally, it is very difficult to access the low temper¬ 
ature region where statistical errors become important. 
As a real frequency method, NRCP can avoid the difficul¬ 
ties that arise from the need to carry out analytic contin¬ 
uation. However, the method becomes extremely cum¬ 
bersome for more than one impurity or channel. NRG 
is better suited for low temperature studies. Recently,!^ 
NRG was applied to study degenerate multi-orbital lat¬ 
tice problems, but the non-degenerate case remains un¬ 
explored. EEEI is also a real-frequency method, but one 
considers only a finite number of bath states, so the re¬ 
sulting energy spectrum is discrete, and the broadening 
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procedure for obtaining continuous spectra is not free of 
ambiguities. Moreover, large systems or multi-orbitals 
are not accessible. DMRGP for the single site case has 
some numerical artifacts and its accuracy as an impurity 
solver is not entirely cleaJ^. 

The semi-analytical methods are perturbation the¬ 
ory based solvers that attempt to capture the essential 
physics by constructing an ansatz for the single-particle 
quantities. The ansatz is based on satisfying various lim¬ 
its or conservation laws, and comprises diagrams up to 
a certain order or sums a specific class of processes to 
infinite order. The main advantages of these methods 
are that they are computationally less expensive than 
the numerically exact methods listed above, while also 
yielding real frequency data. However, semi-analytical 
methods are, by definition, approximate and need to be 
benchmarked against exact results to gauge their range of 
validity. For example, although NCA^ gives qualitatively 
correct results for temperatures higher than the Kondo 
temperature, spurious non-analyticity at the Fermi en¬ 
ergy develops at lower temperature^. To recover the 
correct Fermi liquid behavior at low temperatures, one 
needs to consider a larger class of diagram^^Sl xhe FLEX 
approximation is conserving in the Baym-Kadanoff sense, 
but it does not have the correct strong coupling behav¬ 
ior. So when it is employed for the half-filled Hubbard 
model, strong coupli ng physics like the Mott transition 
is not capture d^^b^ l The FLEX[i^ has been extended 
to study degenerate multi-orbital problems but the is¬ 
sues plag uing s ingle-orbital problems remain. The LMA 
is a highljISEH accurate technique, that has been bench- 
marked extensiveljffil with NRG, but the method has not 
been used to study lattice problems except the periodic 
Anderson modeP^. Moreover, extensions to symmetry 
broken phases or multiple orbital problems remain to be 
carried out. 

The IPT is a simple, second order perturbation the¬ 
ory bas ed me thod and it has been used widely to solve 
impurit} f^^b9 | and lattice problem^^ at zero as well as at 
finite temperature. In the IPT, a self-energy ansatz is 
constructed that interpolates between known limits (i.e., 
weak coupling, atomic and high frequency limits) which 
is why it is also called an interpolative approach. It is 
clear that even the single-orbital IPT is not free of ambi¬ 
guities so different constraints or limits to construct the 
ansatz yield different results. Hence, an IPT for multi¬ 
orbital problems has been ‘synthesized’ in many different 
ways by various group^^, and we discuss these variations 
next. 

The IPT ansatz for the self-energy E(a;) is based on a 
rational or continued fraction expansion of a specific sub¬ 
set of diagrams, and consists of a small number of free 
variables that are fixed by satisfying various limits, such 
as atomic and high frequency limits and conservation 
laws such as the Luttinger’s theorem. Such an inter pola- 
tive approach was first initiated by Martin RodercP^ED 
for the single impurity and periodic Anderson models. 
The approach used the second order self-energy as a 


building block and the pseudo-chemical potential fiQ, was 
fixed by assuming that the occupation no of the non¬ 
interacting part of the Anderson impurity problem is 
equal to the lattice occupation n. Soon after the devel¬ 
opment of DMFT, the single band Hubbard model was 
studied by Muller-HartmanrP^ using self-consistent per¬ 
turbation theory E[G]. The self-consistent perturbation 
method is able to produce a coherence peak in the single 
particle spectral function; however, it fails to reproduce 
the high energy Hubbard bands. For the single impu¬ 
rity Anderson model (SIAM), Yosida and YamadcP^^ 
demonstrated that perturbation theory in U is quite 
well behaved for the symmetric case when expanding 
around the Hartree-Fock solutio n. Based on these find¬ 
ings, Georges and Kotliail^Sl^si introduced an impurity 
solver called iterative perturbation theory (IPT) to solve 
the single band Hubbard model within DMFT which is 
based on mean-field Gq and is able to capture coher¬ 
ent and incoherent features of single particle spectrum 
quite successfully. This is one of the biggest advantages 
of theories based on Go rather than the fully dressed 
G. Another advantage they offer is that, these theories 
naturally avoid any form of two particle divergence and 
are therefore able to provide a reliable description also 
“beyond” the perturbative regim^^. This is very impor¬ 
tant, because the twoparticle divergencies, which might 
induce ambiguities in the numerical determination of the 
DMFT self-energy (E[G]) within standard perturbation 
theory schemes, occur in a rather large portion of the 
phase diagram, includiiM the metallic regime much be¬ 
fore the Mott transition^. 

Subsequently, Kajueter and Kotliail^Sl^ proposed a 
modification to the IPT called modified iterative per¬ 
turbation theory (MIPT). In addition to the usual con¬ 
straints of IPT, the MIPT constrains the zero frequency 
behavior of the self energy by adding a pseudo chem¬ 
ical potential /ig to the Hartree corrected bath propa¬ 
gators. This pseudo-chemical potential, /rg, can be ob¬ 
tained in different ways so there is an ambiguity in the 
method. Kajueter!^ fixed this free parameter by satisfy¬ 
ing the Friedel’s sum rule (equivalently Luttinger theo¬ 
rem), hence his method is called IPT-L. The Luttinger 
theorem and Friedel’s sum rule are valid only at zero 
temperature, hence for finite temperature calculations, 
Kajuetei!^ used the same /ig that was obtained at zero 
temperature. 

To study spontaneous magnetism in the si ngle b and 
Hubbard model, Potthoff, Wegner and Noltin^^oiMJ j]gg_ 
proved MIPT further by taking into account the spectral 
moments up to third order and instead of fixing /ig by 
using Luttinger theorem, they fixed it by the n = ng 
constraint. This method may be called IPT-ng. They 
also considered the simpler option, where lattice chemi¬ 
cal potential /i is equal to the pseudo chemical potential 
/ig. This is called IPT-/ig and they bench-marked IPT-L 
and IPT-ng with IPT-/ig. Recently, Arsenault, Semon, 
and TremblajEH bench-marked IPT-ng with CTQMG and 
found the pathology in IPT-ng that, in the strong cou- 
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pling regime, the method does not recover a Fermi-liquid 
for filling close to n = 1. They suggested a new method 
(IPT-D) to fix the through a double occupancy con¬ 
straint. The range of schemes originating from the inher¬ 
ent ambiguities at the single-orbital level give an idea of 
the far larger range of approximations that can be built at 
the multi-orbital level. These schemes will be described 
next. 

KajueteJ^ extended his single orbital perturbative 
scheme to the degenerate multi-orbital case. He used 
the coherent potential approximation (CPA) to calculate 
higher order correlation functions in the self energy. He 
showed, by benchmarking against ED, that the scheme 
provides reasonable results only if the total particle den¬ 
sity per site is less than one. For fillings greater than one, 
his scheme produced a false double peaked structure at 
the Fermi level instead of a single resonance. The reason 
for such a spurious structure is that the high frequency 
tails in the continued fraction expansion can be system¬ 
atically improved by considering poles involving higher- 
order correlations functions in the self-energy, but this in 
turn seriously degrades the low frequency behavior when 
the Luttinger’s theorem is attempted to be satisfied. To 
study quantum transport in mesoscopi c sy stems such as 
multi-level quantum dots, Yeyati et al.*^ introduced an 
interpolative scheme based on IPT-riQ. LiebsctP^ ap¬ 
plied an extension of IPT to study the orbital selec¬ 
tive Mott-transition, using which he showed that inter¬ 
orbital Coulomb interactions gives rise to a single first- 
order transition rather than a sequence of orbital selec¬ 
tive transitions. In Liebsch’s extension of IPT for the 
multi-orbital case, he chose the self energy to be the com¬ 
bination of Hartree term and second order pair-bubble 
diagram with interaction vertices between electrons in 
different orbitals on the impurity. Laad et al.^^ con¬ 
structed an interpolative scheme for multi-orbitals that 
was used extensively to study real materials through the 
LDA-I-DMFT framework. In a similar context, Fujiwara 
et al.l^ developed an interpolative approach for degener¬ 
ate multi-orbitals. The novelty of their method was that 
they used ligand field theory in the atomic limit to find 
the higher-order correlation functions. 

Although there exist a large range of schemes for ex¬ 
tending IPT to the multi-orbital case, extensive bench¬ 
marking of any single method has not been carried out. 
Recently Savrasov et al.l^, and Oudovenko et al.l^ de¬ 
veloped an interpolative approach for degenerate multi¬ 
orbitals based on a simple rational form of the self-energy, 
where the unknown coefficients in the self-energy are de¬ 
termined using slave boson mean-field and Hubbard I 
approximations. In their Hirsch-Fye-QMC work on the 
SU(4) Hubbard model, they have observed a good agree¬ 
ment in the particle-hole asymmetric cases. 

In the present work, we build upon the previous knowl¬ 
edge to develop an interpolative scheme for solving a 
general multi-orbital quantum impurity problem. Our 
scheme is also based on the second-order self-energy 
as a building block and we use the generic name for 


the method as simply multi-orbital iterative perturba¬ 
tion theory (MO-IPT). Our method has a single pseudo¬ 
chemical potential that is found by satisfying the Lut¬ 
tinger’s theorem. We impose the correct high frequency 
and atomic limits to get the unknown coefficients in the 
self-energy ansatz. In the single orbital case, we find that 
MO-IPT recovers the usual MIPT self energy expression 
and for the degenerate multi-orbital case, our MO-IPT 
self-energy expression reduces to that of Kajuetei!^. The 
main novelty lies in handling the high frequency poles in 
a systematic way. The method is general enough that it 
can be applied to study symmetry broken phases, Hund’s 
coupling (density-density type) and crystal field effects. 

Since MO-IPT is a semi analytical method it needs 
to be bench-marked. Subsequent to the description of 
the method, we embark upon an extensive benchmark¬ 
ing of MO-IPT with numerically exact, hybridization ex¬ 
pansion continuous time quantum Monte Carlo method 
(S-CTQMC)P^ as implemented in the ALP^^ libraries 
and our own implementation of interaction expansion 
CTQMC (W-CTQMC). Our main conclusion is that the 
MO-IPT method works very well when used away from 
integer-fillings, even at reasonably strong coupling. Us¬ 
ing MO-IPT, we have addressed issues disputed in the 
current literature of doped Mott insulator^ and cova¬ 
lent band insulator^Sl jn the multi-orbital degenerate 
case, the method proposed by Kajueteil^ shows spuri¬ 
ous features which restricted its applicability to fillings 
smaller than one. However, our approach circumvents 
all the above issues and moreover captures the filling de¬ 
pendent effect of the Hund’s couplinj^ in the low en¬ 
ergy scale. Our study of crystal field effects (the non¬ 
degenerate case) is the first attempt to extend the ansatz 
beyond the degenerate case. In addition, we have shown 
that our method produces a good default model for the 
analytic continuation of CTQMC data using the max¬ 
imum entropy methocP^. We have also integrated the 
MO-IPT with material-specific, density functional the¬ 
ory based calculations, and thus tested it for a proto¬ 
typical example of strongly correlated electronic system, 
SrVOa. A rather good agreement is obtained when the 
MO-IPT photo-emission spectra (PES) is compared with 
experiments. 

We have organized the paper as follows. In section 
H, we outline the formalism for MO-IPT. In section HI, 
we discuss results when the MO-IPT is applied in the 
DMFT context for lattice problems. In section IV, we 
present our conclusions along with future directions and 
possible improvements. 


II. Model and Formalism 


The multi-orbital Hubbard model for density-density 
type interactions and for cubic environment in standard 
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second quantization notation is given by 
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( 1 ) 


-ff(k)) ^ is the hybridization matrix or equivalently the 
self-consistently determined bath and Simp(w^) is the 
Tj^^mpurity self-energy obtained by solving the impurity 
problem. The set of equations is closed by noting 
that, within DMFT, the lattice self-energy is momentum- 
cr. independent and is the same as the impurity self-energy, 
i.e S(k, w'*') = Simp(w^). The local Green’s function ob¬ 
tained in Eq. (§ is used for defining a new hybridization 
as 


where creates the electron at lattice site i, in orbital 
a with spin a and Cj^a annihilates the electron at site j, 
in orbital /3 with spin a. We are mainly interested in the 
local single particle electron dynamics, which is given by 
the momentum sum of the lattice Green’s function 

Gioc{uj+) = Y -^^-• (2) 

^(a;++M)I-i?(k)-S(k,cc+) 

Where w"*" = uj+'ir] and r] is the convergence factor. Here 
iJ(k) comprises intra-unit-cell hybridization and inter¬ 
unit-cell hopping, namely 

i?(k) 

— Hintra + i?(k) 

inter ( 3 ) 

where (Hintra) = (4) 

\ / Q;/3 

and =e(k)„^, (5) 

where e^o, are orbital energies, T“j^ are intra-unit cell hy¬ 
bridization matrix elements, and e(k)Q ,/3 is the dispersion 
of the lattice, that depends on its geometry. For exam¬ 
ple, in the case of a simple cubic lattice, e(k)Q,^ assumes 
the form, — 2f“^ (cos kx + cos ky -\- cos kz). 

Within DMFT, one can map the multi-orbital Hub¬ 
bard model on to an auxiliary impurity problem with 
a self consistently determined bath. The Hamiltonian of 
the corresponding single impurity multi-orbital Anderson 
model, is expressed in standard notation as: 

T~Limp = ~ ^ E! ^ka -|- 

ot k^OL 

+ E + pUanp (6) 


Here a and /3 are impurity orbital indices including spin. 
The first term in the above equation represents the or¬ 
bital energy; the second term is the hybridization be¬ 
tween the impurity and the host conduction electrons, 
the third term represents the host kinetic energy and 
the final term is the local Goulomb repulsion between 
electrons at the impurity. The corresponding impurity 
Green’s function is given by. 


Gimp 


(w+ + n)l-e- A(w+) - Sjmp(w+) 


( 7 ) 


where (e)a /3 = Ca^a/?- Here A(a;+) = |Vfcap(a;+I- 


A(w+) = (a;+ +fi)l-e- t^mpY) - G-^Y ■ (8) 


Obtaining the self-energy however is the most challeng¬ 
ing step, and we employ multi-orbital iterated perturba¬ 
tion theory to solve the multi-orbital Anderson model. 
The starting point, as usual, is an ansatz for the impu¬ 
rity self-energy, given bji^ 



^aj3 


I 'y ^ Uajin-f) + 

\77^“ 




( 9 ) 


The self-energy is thus restricted to being diagonal in 
the orbital basis. In the above ansatz, the first term is 
simply the Hartree energy and the second term contains 
the second order pair-bubble diagram of matrix size 
NxN, where N is the number of orbitals. The second 
order pair-bubble diagram on the real frequency axis is 
given by 


.( 2 ) 


=Ulp J J J deide2de3pa{ei)pfi{€2)pfi{e3) 

nF(-ei)n_F(e2)nF(-e3) -|- nF(ei)nF(-e2)nF(e3) 


OJ 


+ _ 


ei -I- £2 - £3 


( 10 ) 


where pa = and Qaa is the Hartree corrected 

bath propagator, which is obtained from a Dyson like 
equation, and is given by 

(GroES + e-(M-A^o)l) . (11) 

The pseudo chemical potential, po, is found at T = 0 by 
satisfying the Luttinger’s theorem. 

At finite temperature, an ambiguity exists in the deter¬ 
mination of the pseudo-chemical potential. We choose to 
use the po determined at zero temperature for all finite 
temperatures. The chemical potential, p, is found by fix¬ 
ing the total occupancy from the local Green’s function. 
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Gioc, to be equal to the desired filling, 

1 

-Im / TrGioc = ritot , (13) 

^ J — oo 

where the trace is over the spin and orbital indices. The 
unknown coefficients Aa , from Eq. 0 are obtained in 
the standard way by satisfying the high frequency limit 
and the atomic limit respectively. The detailed procedure 
to derive A^Ba and their expressions are discussed in 
the Appendix. These coefficients contain higher order 
correlation functions. The order of the correlation func¬ 
tions depends on number of poles in the self energy. For 
example a pole of order n involves (n-fl)*^ order correla¬ 
tion functions. For a two pole ansatz Aa and Ba involve 
two and three particle correlation functions. We calculate 
the two particle correlation functiorP^l using the equation 
of motion method to obtaiiJ^ 


Umni'innirim') J (a;)Im (w)G'm(w)] . 

m'^m 

(14) 

This single equation is not sufficient to find all the two- 
particle correlators. Hence as an approximation, we use 
the following: 


/da;nF(a;)Im(Sm(a;)Gm(w)) 

=--1)—^ 

We calculate the three particle correlation function en¬ 
countered in Bq, approximately by decoupling it in terms 
of two and single particle correlation functions. In this 
work, we have ignored the three particle correlation func¬ 
tion. 


III. Results and Discussion 


The formalism developed in the previous section is ap¬ 
plied to a wide variety of correlated systems. We begin 
with a discussion of the well studied paramagnetic Mott 
transition in the half-filled single-band Hubbard model. 
Then we examine the doped Mott insulator. The covalent 
insulator is considered next, followed by the two-orbital 
Hubbard model. For the latter. We investigate the effect 
of filling, Hund’s coupling and crystal field splitting. Fi¬ 
nally, we move to real material calculations considering 
specifically the case of SrVOs. As mentioned earlier we 
bench-mark our results with those from numerically ex¬ 
act CTQMCP methods. The CTQMC formalism yields 
results on the Matsubara frequency axis so to get the 
real frequency data, analytical continuation is required. 
We avoid analytic continuation by transforming the real 
frequency data obtained from MO-IPT to imaginary fre¬ 


quencies using the following spectral representations: 


f AG{u})duj 

G{iuin) = / ^-, (16) 

J lUJn- ^ 

and 

f AY:{u})duJ 

S(iw„) = / ^-, (17) 

J lUJn-OJ 

where Ag(w) = —ImG(w)/7r and As(w) = —ImS(a;)/7r. 
In order to quantify the efficiency of the method, the 
imaginary part of the self energy needs to be bench- 
marked rather than the Green’s function. This is because 
the former is far more sensitive than the latter and more¬ 
over, the low energy scale of the system depends on the 
imaginary part of the self energy. 


A. Single band Hubbard model: Half-fllled case 


The Hamiltonian for the single band Hubbard model 
is given by 

T~L = 'y ( + h.c) + y ( eio-nio- + y ( —niariia . 

ij(7 i<7 i<j 

(18) 

We study the above model within DMFT for a semi¬ 
elliptical density of states, given by 



Here W is the full-band width. In our calculations, we 
choose the energy unit to be ^ = 1. 

The half-filled Hubbard model exhibits an interaction- 
driven metal-insulator Mott transition at a critical Uc- 
Terletska et al.l^ found that the critical exponents and 
scaling functions obtained by IPT are identical to those 
from CTQMC. Here, we revisit this case and benchmark 
the quasiparticle weight, double occupancy, spectra and 
imaginary part of the self-energy. The MO-IPT method 
reduces to the second order perturbation theory in terms 
of Hartree-corrected propagators. In Figure[^a) we com¬ 
pare the quasi-particle weight Z obtained from different 
impurity solvers and several values of the Coulomb inter¬ 
action. The values of Z obtained from S-CTQMC match 
well with those from NRGP^ for all values of U/W except 
close to the Mott-transition. This is most likely because 
we have done CTQMC calculations at /I = 64, while NRG 
is at zero temperature. The critical interaction strength, 
^ « 1.35 obtained from both the method^ agrees very 
well. The Z obtained from MO-IPT at ,5 = 64 matches 
quantitatively with CTQMC and NRG in the weak cou¬ 
pling limit and only qualitatively in the proximity of the 
transition. On the other hand, the results of the self en¬ 
ergy functional approach (SFA)!^ agree with MO-IPT in 
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FIG. 1. (Color online) (a) Quasi-particle weight Z of the 
single band half-filled Hubbard model obtained with differ¬ 
ent impurity solvers (see text for more details) (b) Double 
occupancy D obtained from MO-IPT and S-CTQMC. 


the strong coupling limit rather than in the weak cou¬ 
pling limit. The MO-IPT yields the critical value of 
^ = 1.42, which is in good agreement with the critical 
value ^ = 1.45 obtained from SFAESI at zero tempera¬ 
ture. The double occupancy obtained from MO-IPT and 
S-CTQMC (shown in panel (b) of Fig. also match, 
except very close to the transition. A detailed compar¬ 
ison of spectra from S-CTQMC and W-CTQMC with 
the same from MO-IPT (transformed to imaginary fre¬ 
quencies) is shown in Figure]^ The left panels show the 
imaginary part of the Green’s function at U/W = 1.0 
(top panel) and U/W = 1.5 (bottom panel), while the 
right panels show the imaginary part of the correspond¬ 
ing self-energies. The excellent agreement between the 
three methods is clearly evident. 


B. Single band Hubbard model: Doped Mott 
insulator case 

The single band Hubbard model has gained a lot 
of interest, because the doped Mott insulator regime 
is believed to capture the essential physics of high Tc 
superconductorsSB. This regime is, in reality, highly 
complex, because many different factors such as prox¬ 
imity to the antiferromagnetic Mott insulator, disor¬ 
der, d-wave superconducting fluctuations and pseudogap 
physics have to be treated on an equal footing. Hence, 
investigations of the doped Mott insulator in all its glory 


FIG. 2. (color online) Gomparison of the imaginary part 
of Matsubara Green’s function (left panels) and self energy 
(right panels) obtained from MO-IPT, S-CTQMG and W- 
GTQMCPfor U/W = 1.0 (top panels) and U/W = 1.5 (bot¬ 
tom panels) at /3 = 64. 

represents one of the toughest challenges in condensed 
matter. Here, we take a simplistic approach to the prob¬ 
lem, and investigate the performance of MO-IPT in the 
paramagnetic doped Mott insulator in infinite dimen¬ 
sions. Our MO-IPT reduces basically to the IPT-L in 
this regime. 

A comparison of quasi-particle weight at U/W = 1.5 
obtained from MO-IPT and S-CTQMC as a function of 
filling (Fig. ^ yields, surprisingly, an excellent agree¬ 
ment. We observe that as we decrease the filling (from 
1) for a given U/W, the Mott insulator turns into a 
strongly correlated metal and finally ends up as a sim¬ 
ple metal. In the strong coupling limit, for filling close 
to n = 1, the IPT-no method gives an insulating solu¬ 
tion, while the IPT-L correctly predicts a metal in agree¬ 
ment with exact methods. Kajueter and Kotliail^ have 
benchmarked the real-frequency spectral functions ob¬ 
tained from IPT-L with exact diagonalization calcula¬ 
tions and had found good agreement. We find that the 
imaginary part of the Green’s function and self-energy 
obtained from IPT, when transformed to the Matsubara 
frequency axis using Eqs. I®, (0 are almost identical 
to those obtained from the strong coupling and weak- 
coupling variants of CTQMC (see Fig. [^. The slope of 
the ImE(iw„) as a;„ —?► 0 is 1 — \/Z, and the good agree¬ 
ment of Z shown in Fig. is simply a reflection of the 
detailed agreement for all frequencies. Such an excellent 
agreement is truly surprising because IPT is a perturba¬ 
tive method by construction and the strongly correlated. 
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FIG. 3. (color online) Quasi-particle weight obtained from 
MO-IPT (or IPT-L) is compared to the same obtained from 
CTQMC for the paramagnetic doped Mott-insulator as a 
function of filling with VjW =1.5 and (5 = 64. 



FIG. 4. (color online) Doped Mott insulator: Gomparison of 
imaginary part of Matsubara Green’s function and self energy 
obtained from MO-IPT, W-GTQMC and S-CTQMG for U/W 
= 1.5 at different fillings and /I = 64. 


doped Mott insulator regime should not, in general, be 
amenable to perturbative methods. 


C. Covalent Insulator: 


The discovery of topological insulator^^ has led to a 
renewed interest in the role of e-e correlations in band 
insulators (BI)pl. The prime examples of such materials 
would be FeSPl and FeSbpl, since experimental mea¬ 
surements indicate a small optical gap and large ther¬ 
mopower (at low T). Increasing temperature leads to 
closing of the gap, and concomitantly a insulator-metal 
crossover in the resistivity. Such large scale spectral 
weight transfers are highly indicative of strong corre¬ 
lations. Specific heat measurements also seem to vali¬ 
date this observation. The band gap in these systems 
is a simple consequence of the structure of the hopping 
matrix and not of completely filled electronic sh ells^. 
Hence these materials are called covalent insulator^SMl 
A Hamiltonian that describes the covalent insulator is 
given bjna 

«=E (“U 4.,) A(k) 

k(T ^ 

( 20 ) 

where a = a and b are two sub-lattices with semi-elliptic 
bands and having dispersion Ck and -Ck respectively. The 
two sub-lattices are coupled by a k—independent hy¬ 
bridization V. While the unit of energy is chosen to be W 
= 2 throughout, for this subsection W=4 has been chosen 
in order to benchmark with earlier result^^. This is the 
first two-band model we have studied in this work, since 
the previous cases pertained to the single-band Hubbard 
model. Hence this will be the first real test of the ‘multi- 
orbital’ part of MO-IPT. Since this is still the half-filled 
case, Luttinger’s theorem does not have to be satisfied 
explicitly. The = 1 and Hq, = 0 for all orbitals. 
Thus, the MO-IPT used for the covalent insulator case 
is equivalent to that employed by LiebsclP^ for studying 
the Mott transition in the two-band Hubbard model. 

Before presenting the numerical results, we will present 
a simple analytical argument about the absence of any 
intermittent metallic phase in the CBI at zero tempera¬ 
ture unlike the interaction-induced metallic phase in the 
Ionic Hubbard modeP^. The structure of Ho-(k) in the 
CBI is given by, 


E 


u, 


O'ca'^ia'Ha 7 


H„(K) 




( 21 ) 


and the corresponding impurity Greens function for sub¬ 
lattice ’a’ is given by, 



, (ba(uJ'^,e)po(e) 
^Ca^(w+,e)Cba(w+,e) - F2 


( 22 ) 


























where 


Coo-(w'', e) = uj + iT] + fi-e - Eao-(w+), 
e) = uj + ir] + n + e- , 


and 77 —?► 0'*'. In the half-filling case, the Hamiltonian has 
mirror type symmetry between sublattices, which reflects 
in the impurity Green’s function and self-energy in the 
following way, 

= - [GtA-^^)]* , (23) 

= U - ■ (24) 


By using above self-energy symmetry relation, it is easy 
to see that 


Cao'(a7 , e) — [Cbo’( ^ J* 


(25) 


then Eq. (22) can be written as. 


Gaa{0O+) = J de 




Ca.(cc+,e)G.(-w+,e) + 4^^ ’ 


(26) 


With the assumption that in the band insulator and 
metallic phases, a Fermi-liquid expansion of self-energy 

holds, namely that E(a;) Ti{Q)+uj{1 — 1/Z)-\-0{ijj'^). 
Then, the value of imaginary part of self-energy at zero 
frequency (and zero temperature) is ImEao'(O) = 0, and 
the corresponding density of states (DOS) Daa{Q) = 
-ilmGaa(O) is given by, 


D (0) = f __ /2J) 

7 ^2 + [Re(C„,( 0 ,e ))]2 + E 2 ’ ^ 

where Re(Cao'(0, e))=[ 7 i — e —ReEacr(O)]. In the limit 77 —>• 

0 +, 

Daa(0) = J depoie)S + [p - ReEa< 7 ( 0 ) - e)^) . 

Since the argument of the Dirac delta function is posi¬ 
tive definite for any V ^ 0, the density of states at the 
chemical potential, Daa{Q) will necessarily vanish for any 
R 7 ^ 0, and hence the system will be gapped. Thus for 
the CBI, interactions do not close the gap, no matter how 
strong they are. This implies a clear absence of metallic- 
ity in these insulators. 



FIG. 5. (color online) Covalent insulator: (a) Quasi-particle 
weight Z as a function of U/W obtained from MO-IPT (black 
circles) and CTQMC (red squares) for /3 = 60 and V=0.5.(b) 
Double occupancy as a function of [7/IT obtained from MO- 
IPT and S-CTQMC. (c) Charge gap as a function of U/W 
obtained from MO-IPT at T=0. 


gap with decreasing temperature. Precisely this behavior 
is seen in the real frequency spectra (left panels. Fig. |^, 
which arises from the spectral weight transfer in the self¬ 
energy as a function of temperature. The high reliabil¬ 
ity of these spectra and self-energies computed through 
MO-IPT is apparent in the excellent agreement with the 
same obtained through strong coupling CTQMC (on the 
Matsubara axis, in Fig. [^. The crossover of the band- 
insulator to Mott insulator is also visible in the increas¬ 
ing (negative) slope of the imaginary part of self-energy 
with increasing U/W. Similar kind of correlation induced 
transitions have also observed bet ween t opologically triv¬ 
ial and non-trivial band insulator^SlIlll_ 


D. Two orbital Hubbard model: 


The quasi-particle weights (Fig.j^a)) and double occu¬ 
pancy (Fig.j^b)) of sublattice ’a’ obtained from MO-IPT 
and S-CTQMC (shown as black circles and red squares 
respectively) are in close agreement except in the prox¬ 
imity of the transition of the correlated band insulator 
to a Mott insulator. Unlike the ionic Hubbard model 
cas^^, we do not see any intervening metallic phase be¬ 
tween the correlated band insulator and the Mott insu¬ 
lator. This is also consistent with the S-CTQMC and 
analytical results. At high temperatures, the correlated 
band insulator should be gapless, and must develop the 


Encouraged by the excellent benchmarking of MO- 
IPT with CTQMC for the two-band covalent insula¬ 
tor system, we now move on to the two-orbital Hub¬ 
bard model^^lSil Hamiltonian, in standard nota¬ 

tion, for a cubic environment and for unbroken spin sym¬ 
metry, is described in Eq. 0 - Throughout the paper, 
we have considered local interactions of density-density 
type which are obtained by neglecting spin flip and pair¬ 
hopping terms that must be present for a rotationally 
invariant Hund’s coupling. The hopping is taken to be 
diagonal in orbital indices for simplicity. 
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FIG. 6. Covalent insulator: Spectral functions (left panels) 
and imaginary part of self energy of sublattice a (right panels) 
from MO-IPT at UjW = 1.25 and V=0.5 for a range of /3 = 
1/T values (increasing T from top to bottom). 



FIG. 7. (color online) Covalent insulator: Comparison of the 
imaginary part of Matsubara (a) Green’s function and (b) 
self-energy of sublattice ’a’ obtained from MO-IPT (black) 
and S-CTQMC (red) for various UjW values and /3=60. 


(a) Half-filling: J = 0 

We begin by considering the half-filled case (total occu¬ 
pancy is two) with J = 0. The Hamiltonian (Eq. 0 ) has 
SU(4) symmetry in this situation. We have employed a 
semi-elliptic non-interacting density of states of full-band 
width W = 2 for the MO-IPT-DMFT calculations. 



FIG. 8. (color online) (a) Two-orbital SU(4) symmetric Hub¬ 
bard model at half-filling: Quasi particle weight as a function 
of U/W obtained from different impurity solvers(SFA,MO- 
IPT,FLEX,S-CTQMC at /3 = 64 and ED at T=0.0).(b) Dou¬ 
ble occupancy obtained from MO-IPT (black circles) and hy¬ 
bridization expansion CTQMC (red squares) for /3=64. 

In Fig. [^, we plot the quasi-particle weight {Z) ob¬ 
tained from different impurity solvers for the particle- 
hole symmetric case. The results from strong cou¬ 
pling CTQMC, EElSll and SFAI^SI, including the critical 
Uc, where the system transitions from metal to Mott- 
insulator, are in good agreement. The critical value Uc 
obtained in the multi-orbital case is greater than the 
value obtained in the single band case. The Mott tran¬ 
sition is absent in the FLEX result!^. The MO-IPT is 
seen to underestimate the quasiparticle weight as com¬ 
pared to the other methods (except MO-IPT<nn>=0; 
see below). However, the critical Uc agrees reasonably 
well with that from hybridization expansion CTQMC. 
The green diamonds are from a variant of MO-IPT (used 
e.g. by Fujiwara et al.l^ where the two-particle correla¬ 
tion function is simply decoupled into two single-particle 
terms {{nani 3 )={na){nij)). The neglect of two particle 
correlations leads to a much worse comparison than MO- 
IPT. In contrast to the not-so-good agreement with exact 
methods for the quasiparticle weight, the average double 
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occupancy obtained from MO-IPT shows excellent agree¬ 
ment with CTQMC (see Fig. |^). Since the total energy 
of the system depends on single particle and two parti¬ 
cle correlation functions, we expect that thermodynamic 
quantities like total energy or specific heat computed 
through MO-IPT might be reliable. One more impor¬ 
tant observation is that the double occupancy remains 
finite and almost constant even beyond the Mott transi¬ 
tion, unlike the the single band case. We also compare 
the single-particle Green’s function and self-energy on 
the Matsubara frequency axis (Fig. [^. At high frequen¬ 
cies, the agreement between MO-IPT and S-CTQMC is 
seen to be excellent, while the agreement worsens at low 
frequencies, especially with increasing UjW. 



FIG. 9. (color online) Two-orbital SU(4) symmetric Hubbard 
model at half-filling: Imaginary part of Matsubara Green’s 
function (left panels) and self energy (right panels) obtained 
from MO-IPT (red solid lines) and S-GTQMC (black solid 
lines) at b=64. 


(b) Half-filling: Effect of Hund’s coupling (J) 


The interplay of Hund’s coupling, J, and local inter¬ 
action, ?7, has been investigated by several groups. The 
main consensus is that strong corre lation effects can be 
affected significantly through jElEII pgr example, in the 
half-filled case, the C4 for Mott transition is lowered by 
{N — I) J, where N is the number of orbitals, while the 
critical U is enhanced W 3J in the non-half-filled (but in¬ 
tegral occupancy) case^. It is important to know the ex¬ 
tent to which the interplay between J and U is captured 
by the MO-IPT method. In Fig. 10 a) and (b), the quasi- 



FIG. 10. (color online) Two orbital half-hlled Hubbard model, 
finite J. Quasi particle weight dependence on UjW obtained 
from (a) strong coupling GTQMC, (b) ED and (c) MO-IPT 
for various J values. Insets in the panels (a) and (b) show the 
effect of J on Z in the weak coupling regime. 


particle weight Z for different values of J obtained from 
S-CTQMC and EEI^ is shown. Indeed, with increasing 
J, the Uc at which Z —0 decreases sharply, as expected 
from the atomic limit. Also, for each J, the quasiparticle 
weight decreases monotonically with increasing interac¬ 
tion strength. Although the latter trend is qualitatively 
captured by the MO-IPT result (shown in Fig. 10:) for 
larger J, there is a disagreement with the exact results at 
lower J values. The MO-IPT yields a C4 that is a non¬ 
monotonic function of J. The insets of panels a and b in 
Fig. 10 zoom in on the low interaction {UjW < 0.3) part 
of the main panels. Unlike for U/W > 0.3, where increas¬ 
ing J leads to a monotonic reduction of Z, a rise and fall 
of Z is observed for U/W < 0.3. Although such a trend is 
achieved by MO-IPT as well, the non-monotonicity sus¬ 
tains even for larger UjW. A frozen local-moment phase 
is seen in the S-CTQMC calculations for any given J in 
the strong coupling limit, while such a phase is not ob¬ 
served either by EEEsI or MO-IPT calculations. It must 
be mentioned here that the CTQMC calculations employ 
a density-density type Hund’s coupling, while the ED 
employs a fully rotationally invariant J. Although the 
quasiparticle weight dependence on U and J is not accu¬ 
rately captured by MO-IPT, the single-particle dynamics 
on all scales is in qualitative agreement with S-CTQMC 
calculations (as seen in Fig. II). 
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FIG. 11. (color online) Two orbital half-filled Hubbard model, 
finite J. Imaginary part of Matsubara Green’s functions 
(left panels) and self-energy (right panels) obtained from S- 
GTQMG (black) and MO-IPT (red) for different values of J 
and UjW at ,5=64. 


(c) Away from Half-filling: Effect of J 


FIG. 12. (color online) Two-orbital Hubbard model: Effect 
of J away from half-filling {ntot = 1.1). The imaginary part 
of the Matsubara self-energy for various J-values, and fixed 
U/W = 1 as computed within (a) S-GTQMG and (b) MO- 
IPT. Comparison of quasi particle weight obtained from MO- 
IPT (black circles) and CTQMC (red squares) as a function 
of U/W for (c) J = 0.0, (d) J = [7/4 and (e) J = [7/3.5 for 
P = 64; and (f) as a function of J for a fixed UjW = 1.0. 


The MO-IPT method works best away from half-filling, 
which is consistent with the results of comparisons car¬ 
ried out previously by other group^^. In order to illus¬ 
trate this, here we study the two orbital Hubbard model 
for a total occupation of ntot = 1.1. The imaginary part 
of the Matsubara self-energy obtained from S-CTQMC 
matches well with that from MO-IPT (Fig. 12, panels 
(a) and (b)), hence the latter does well in this regime. 
This observation is reinforced by the panels (c)-(e), which 
show a comparison of the quasiparticle weights as a func¬ 
tion of U/W for three values of J, namely J = 0, [//4.0 
and [7/3.5. The results of MO-IPT are seen to agree 
very well with those from CTQMC. For most real ma¬ 
terial calculations, the regime considered in this sub¬ 
section is perhaps the most relevant. Hence, accurate 
results from MO-IPT prove its efficacy for integration 
into first-principles approaches. The Hund’s coupling 
and Coulomb interaction have a synergistic effect at half¬ 
filling, while in the doped case, the reverse occur^^. This 
is shown in panel (f) of Fig. 12 where an increase of Z 


is seen with increasing Hund’s coupling at a fixed inter¬ 
action strength. It is quite instructive to study the real 
frequency spectral functions and self-energies as obtained 
from MO-IPT. These are shown in Fig.j^for various val¬ 
ues of interaction strength and Hund’s coupling, J. In 
the absence of Hund’s coupling, the spectrum (shown in 



FIG. 13. (color online) Two-orbital Hubbard model away 
from half-filling. Real frequency spectral functions (left pan¬ 
els) and minus imaginary part of self energy (right panels) for 
various U/W and J values. 
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the left panels of Fig. 13) exhibits spectral weight trans¬ 
fers characteristic of increasing correlation strength: a 
central resonance that becomes sharper, and Hubbard 
bands that grow in prominence with increasing U/W. 
However, at a fixed UjW^ increasing Hund’s coupling 
leads to a reversal of the aforementioned trend, i.e, a 
broadening of the resonance and a melting of the Hub¬ 
bard band (see e.g. left panel bottom figure). In this 
parameter regime, a previous formulation of the multi¬ 
orbital iterated perturbation theorjEl found a double 
peak structure at the chemical potential. Such a feature 
was shown by the author^^to be spurious by comparison 
to results from exact diagonalization. The reason we do 
not observe such a spurious feature is that we have con¬ 
sidered only two poles in the self-energy, in contrast to 
the formulation of Ref. UHl where they have retained all 
the eight poles (for a two-orbital model). Although our 
ansatz seems like an ad-hoc truncation scheme, the justi¬ 
fication for such a scheme lies in its excellent agreement 
with CTQMC results (shown in Fig. 14) and the absence 
of spurious features. In Fig. the imaginary part of 
Matsubara Green’s functions and self energies obtained 
from MO-IPT are compared with those from CTQMC 
for three values of J at UjW = 1.25 and /3 =64. For all 
values of the Hund’s coupling, an excellent agreement is 
obtained. 



FIG. 14. (color online) Two-orbital degenerate Hubbard 
model away from half-filling (ritot = 1.1). Comparison of the 
imaginary part of the Matsubara Green’s function (left pan¬ 
els) and the self energy (right panels) obtained from MO-IPT 
and S-CTQMC for various values of J at UjW = 1.25. 


E. Two orbital Hubbard model: Crystal field 
splitting and Hund’s coupling 


We now proceed to the case of a non-cteenerate two- 
orbital model with crystal field splitting^ in the pres¬ 
ence of Hund’s coupling. In most materials, the crys¬ 
talline environment lifts the orbital degeneracjl^. For 
example in transition metal oxides, due to crystal field 
effects, the five fold degenerate d-level splits into triply 
degenerate t 2 g and doubly degenerate Bg levels and the 
corresponding energy gap is ^1-2 eV. The degeneracy 
of each of these levels (t 2 g, Bg) is further lifted by dis¬ 
tortions such as the ones occurring in GdFeOa, or aris¬ 
ing through the Jahn-Teller effect or spin-orbit coupling. 
The energy cost for such distortion induced splitting is a 
few meV. Recently, Pavarini et al.^^ studied crystal field 
effects in d^ type perovskites such as SrVOa, CaVOa, 
LaTiOa and YTiOa. It was found that crystal field ef¬ 
fects and cation-covalency (GdFeOa -type distortion) lift 
the orbital degeneracy and reduce the orbital fluctua¬ 
tions. Thus, inve stigating crystal field effects in model 

Hamiltonian^SIlSai jg highly relevant for understanding of 

real materials. 


We have investigated the Hamiltonian in Eq. 0 by 
considering two orbitals with energies ei = 0.0 and 
£2 = —0.2IT, which corresponds to a crystal field split¬ 
ting of 0.2W. The results from MO-IPT, for a fixed to¬ 
tal filling of ritot = 1-1, are compared with those from 
strong coupling CTQMC at the corresponding orbital 
occupancies. In Fig. |15[ we compare the quasi particle 
weights of the two orbitals obtained from MO-IPT with 
that of CTQMC. We observe a better agreement of Z 
for orbital-1 than for orbital-2. This must be expected, 
since orbital-1 is further away from particle-hole symme¬ 
try than orbital-2. The corresponding orbital occupan¬ 
cies as a function of increasing interaction (and hence 
J) are shown in Fig. The deviation between results 


from the two methods increases with increasing U and 
J(= Uj4:), which indicates that MO-IPT is almost exact 
for U/W < 0.5. 


Next, we benchmark the single-particle dynamics in 
the presence of crystal field splitting. In Fig.[l6l we show 
the imaginary part of the Matsubara frequency self ener¬ 
gies obtained from MO-IPT and CTQMC for orbitals-1 
and 2 (left and right panels respectively). The agreement 
between the results is quite evident, this suggesting that 
the MO-IPT should serve as a good method to study in¬ 
teracting, real material systems with finite crystal field 
effects and Hund’s couplings. This is especially true if 
the material in question has a large number of bands, 
which would make it prohibitively expensive to treat with 
CTQMC, while MO-IPT would be able to handle it with 
ease. We now demonstrate the efficacy of MO-IPT when 
applied to a well studied, real material system, namely 
SrVOa. 
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F. Application to real materials: SrVOa 
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FIG. 15. (color online) Crystal field effects. Quasi particle 
weights for (a) orbital-1 and (b) orbital-2, obtained from MO- 
IPT and CTQMC for various U/W values with J = 17/4 at 
/1=64. The insets show the corresponding occupancies. 



FIG. 16. (color online) Crystal field effects. Comparison of 
the imaginary part of the self energy for orbital-1 (left) and 
orbital-2 (right) obtained from MO-IPT and S-CTQMC for 
various values oiUjW and J = 17/3.5. 


Over the past decade or so, the combination of density 
functional theory (DFT) with dynamical mean field the¬ 
ory, such as LDA-|-DMFlt^, has emerged as one of the 
most powerful methods for electronic structure calcula¬ 
tions of strongly correlated electronic systems. Although 
the DFT results contain rich, material specific informa¬ 
tion, being a single particle theory, it works well only for 
weakly correlated systems where the ratio of Coulomb in¬ 
teraction (17) to bandwidth (VF) is small i.e., U/W <C 1. 
If we consider the opposite limit oi UjW ^ 1, we have 
successful methods like the Hubbard-I and Hubbard-III 
approximations or the LDA-I-U method for predicting the 
ground state of the system. But these also have limita¬ 
tions, such as the neglect of dynamical fluctuations in the 
LDA-I-U method. In nature, there are many materials, 
for example transition metal oxides, which lie in between 
these two limits. It has been established in the context 
of model Hamiltonians that the DMFT can handle both 
limits quite efficiently. Hence a natural combination of 
LDA with DMFT is expected to bring predictive capa¬ 
bilities in the theory of strongly correlated electronic sys¬ 
tems. Nevertheless, LDA-I-DMFT is not without its own 
bottlenecks. 

One of the central issues of the LDA-t-DMFT method 
is the correct definition of a correlated subspace. The 
basic idea of a correlated subspace is to make an appro¬ 
priate choice of energy window around the Fermi level 
and fit the band structure to a few-orbital tight-binding 
model. Many techniques have been proposed to con¬ 
struct such a material specific ‘non-interacting’ Hamil¬ 
tonian. The two major techniques for this purpose are 
down-folding^ and projection based Wannier function 
techniqu^^. In general, bands which are crossing the 
Fermi level are considered in the desired energy window 
for Hamiltonian construction. For example in transition 
metal compounds bands having d-orbital character nor¬ 
mally cross the Fermi level. This process becomes simple 
if there is no hybridization in the system and these bands 
with d-orbital character are well separated from other 
bands like the p-bands. As Dang et al.^ pointed out, a 
mixing of these d orbital bands with p orbital bands can 
create several complications. 

After getting the ’non-interacting’ Hamiltonian , one 
can add various types of interactions terms to obtain a 
full material-specific multi-orbital model. The solution 
of such a Hamiltonian is however a major challenge and 
this is where the MO-IPT can be most useful, since it 
scales only algebraically with increasing number of bands, 
while yielding real frequency quantities directly. In con¬ 
trast, impurity solvers like CTQMC and ED scale expo¬ 
nentially with increasing number of orbitals and are very 
expensive, especially for investigations of real materials. 
As a test case, we study SrVOa which is considered a 
prototypical example of a strongly correlated electronic 
system. 
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(a). Computational Details 


We perform our density functional theory (DFT) cal¬ 
culations with linearized augmented plane wave (LAPW) 
based method as implemented in the all-electron package 
WIEN2K!^. The experimentally determined structur^^ 
of cubic SrVOa in a non-magnetic phase was used for 
the calculations (neglecting spin-orbit coupling). The 
product of the plane-wave cut off {K^ax) and the small¬ 
est atomic sphere radius (Rmt) was chosen as Rmt x 
Kmax = 7.0 for controlling the basis set. The radii of 
the muffin-tin spheres were chosen to be 10 — 15% larger 
than the corresponding atomic radii. Thus, the values 
used for Rmt were 2.50 for Sr, 1.89 for V and 1.71 for O. 
With these parameters, charge leakage was absent and 
our DFT results agree well with DFT calculation using 
other basis setPl. We utilize the generalized gradient 
approximation (GGA) of Perdew, Burke and ErnzerhoP^ 
for the exchange and correlation functional. In this calcu¬ 
lation, we consider 512 k-points in the irreducible part of 
the Brillouin zone. After getting the Bloch-eigen states, 
all the necessary inputs for constructing the maximally 
localized Wannier functions (MLWFs) are prepared by 
the WIEN2WANNIER cod^^. Finally, the Hamiltonian 
T-L-dtt is constructed in the maximally localized Wannier 
basis by taking a projection of the three V — t 2 g orbitals 
within the energy window of -1.0 eV to 1.8 eV with the 
standard procedure implemented in WannierQCP^. We 
begin by discussing the DFT results. 


(b). GGA+DMFT: Results and discussion 


Our computed band structure and density of states 
(DOS) are presented in Fig. 17 and Fig. 18 The three 


bands, crossing the Fermi level, are highlighted in cyan, 
violet and grey colors. These bands originate from the 
V — t 2 g states, and are located between -l.leV and 1.5eV. 
The V—Cg states lie at higher energies, between l.leV to 
5.8eV (see the projected density of states in Fig.[T^. The 
band structure agrees well with previous results by Ishida 
et. al.l^ obtained in the LAPW basis. When compared 
with results from the linear muffin-tin orbital (LMTO) 
calculations of Nekrasov et al.^^, the position of P — t 2 g 
bands agrees well but the position of P — states differs 
by about 0.3 eV. This discrepancy is, most likely, due to 
the difference in basis sets used in the two calculations. 
A significant computational simplification results from 
ignoring the hybridization between P — t 2 g and V — Cg 
orbitals, since the low energy correlated subspace com¬ 
prises just the three P — t 2 g orbitals. Thus, the DFT 
results yield a ‘non-interacting’ Hamiltonian 
which in this case is a 3 x 3 matrix for each k. Thus, the 
full DFT + DM FT Hamiltonian is given by 


R = + Hint , 


(28) 



FIG. 17. (color online) Band structure of SrVOs obtained 
from DFT. 



FIG. 18. (color online) The projected density of states (DOS) 
of SrVOa as calculated by GGA (LAPW). 


where Hint is the interaction term is given by 

Hint — R ^ ^ T ^ ^ {U • 

i,a. ioL^^ 

(29) 
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FIG. 19. (color online) Comparison of spectral function of 
SrVOs obtained from different methods for U = 3.44 eV and 
J = 0.46 eV (see text for details). 


In the above expression, i stands for V sites and a is 
the t 2 g orbital index with spin a. U, U'{= U — 2J) 
and U' — J{= U — 3J) are the local, intra orbital and 
inter orbital Coulomb repulsion respectively and J is 
the Hund’s exchange. The local, non-interacting lattice 
Green’s function, in the orbital basis, (G'o(w)), can be 
obtained from the DFT calculated HdftO^) by the fol¬ 
lowing equation as 


Go(w) = 


~ — Hdc 


-1 


(w'*" -I- ^)I — A(w) 


(30) 

(31) 


where fi is the chemical potential and A(a;) is the hy¬ 
bridization. In the DFT approach electronic correlations 
are partially entered through the LDA/GGA exchange- 
correlation potential. This part of the interaction (Huc) 
has to be subtracted in the LDA-l-DMFT approach to 
avoid double-counting. This is not an important issue 
when the low energy effective Hamiltonian contains only 
the d-manifold because we can absorb it into the chemi¬ 
cal potential. However it is an important issue when the 
low energy effective Hamiltonian contains 0-2p orbitals 
also. Various schemes for finding the double-counting 
correction Hjjc exist, each with a different physical mo¬ 
tivation. Details about such schemes may be found in the 
work by Lechermann et al.l^ and Nicolaus ParragtP2lHll. 


In general we can construct the modified host Green’s 
function for the orbital as 


Ga = 


Gq ^ e + Hoc ~ ~ Mo) I 


(32) 


We find the pseudo-chemical potential using the same 
procedure as in the model calculations. The self-energy 

can be found, e.g. through the MO-IPT method out- 

(2) 

lined in Section H The second-order self-energy in 
Eq. (§ is a functional of the modihed host Green’s func¬ 
tions, The full local Green’s function for the lat¬ 


tice Hamiltonian (Eq. (28l) is given by 


G = 


k 


[(a;+ + - HdftO^) - Hoc - E(w) 


n -1 


(33) 

The above Green’s function may be used to obtain a new 
host Green’s function through the Dyson’s equation; 


(/(w) — G -|- E -|- Hdc + £ ~ (m ~ Mo) I 


(34) 


In general, the chemical potential, /i is found by hxing 
the total occupancy from the full Green’s function, G to 
be equal to the value found from DFT, 


Im / TrG = nPJ^ . 


(35) 


where the trace is over spin and orbital indices. 

Thus the full solution of the problem proceeds as fol¬ 
lows. Given the Hdft{^), we guess an initial self-energy, 
as well as the /x and ^o; and use these to find the local 
and the host Green’s functions through Eqs. ( [M| and 
(34). The host Green’s functions are then used to find 


the self-energy, E and Eqs. (33) and (35) are used to find 


the chemical potential. For a fixed /xq, these equations are 
then iterated, until the self-energy converges. With the 
chosen pseudo-chemical potential, the Luttinger’s inte¬ 
gral, Eq. (12) is computed using the converged self-energy 


and local Green’s functions. If the Luttinger’s theorem 
is satisfied within a numerical tolerance, the solution is 
considered to be obtained, else the is tuned, and the 
DMFT equations are iterated, until the Luttinger’s the¬ 
orem is satisfied. 

The DFT predicted occupancy per spin on the three 
correlated V-t 2 g orbitals in SrVOa is 0.166, which implies 
SrVOa is a d^ system. For the DMFT calculations, we 
employ interaction parameters U = 3.44 eV and J = 0.46 
eV, that were obtained by Taranto et al.^^ through the 
random phase approximation (RPA). For SrVOa, we have 
not introduced explicit double counting correction be¬ 
cause we choose the correlated subspace that is identi¬ 
cal with the set of Wannier bands. We absorb the dou¬ 
ble counting correction and orbital energies in the lattice 
chemical potential, which we find by using Eq. (35). 


Our computed GGA-I-DMFT spectrum for SrVOa is 
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shown in Fig. 19 and compared with results obtained 
from other impurity solvers. The GGA result (shown in 
blue) has no signatures of correlation, while each of the 
DMFT calculations exhibit a three peak structure. The 
CTQMG results from GW+DMFT (black) agree quali¬ 
tatively with those from LDA-I-DMFT. However, the de¬ 
tails do differ. Namely, the positions and weights of the 
resonance at the Fermi level and of the Hubbard bands 
differ to a significant extent. This difference, naturally, 
can be attributed to the different starting points, namely 
GW vs LDA, of the CTQMG calculations. Results from 
the MO-IPT solver agree with those from CTQMG in the 
neighborhood of the chemical potential as well as in the 
proximity of the lower Hubbard band. The upper Hub¬ 
bard band is clearly in disagreement with the CTQMG 
results. 



-3 -2.5 -2 -1.5 -1 -0.5 0 

CO(ev) 


FIG. 20. (color online) Comparison of photo emission 
spectra obtained from different methods GW-fDMFT®^, 
GGA-tDMFT(MO-IPT), LDA-tDMFT(CTQMC)P and 
experimenfP^. 

As a final benchmark of the GGA-I-DMFT(MO-IPT) 
calculation, we compare our result with the exper¬ 
imentally measured photo emission spectrum (PES) 
which is shown in Fig. A Hubbard satellite at 

~ —1.5 eV is seen in the experimental PES spec¬ 
trum. Our GGA-I-DMFT(MO-IPT) calculation predicts 
the Hubbard satellite at -1.25 eV. Results from other 
approaches, namely LDA, LDA-fDMFT(CTQMC) and 
GW-I-DMFT(CTQMC) are also reproduced. Surpris¬ 
ingly, the closest match with the experiment is achieved 
by the GGA-I-DMFT(MO-IPT) in terms of the position 
and width of the resonance at the Eermi level and of the 
lower Hubbard band. Thus, we infer that the MO-IPT 
method outlined in this work may be used as an efficient 


tool to study the electronic structure of real material sys¬ 
tems. 


IV. Conclusions 

The development of iterated perturbation theory 
as an impurity solver for single band models and for 
multi-band models dates back to almost two decades. 
Although a few comparisons with numerically exact 
methods have been made, being a perturbative ap¬ 
proach, the method has suffered from reliability issues, 
especially for multi-orbital systems. Nevertheless, sev¬ 
eral multi-orbital extensions of IPT have been proposed 
and used to investigate model Hamiltonians and even 
real material systems. In this work, we have outlined a 
multi-orbital extension of IPT, and benchmarked it ex¬ 
tensively against continuous time quantum Monte Carlo 
results. Our work is the first systematic study of the 
multi-orbital Hubbard model using MO-IPT as a solver 
and varying parameters such as filling, Hund’s coupling, 
and Coulomb repulsion, as well as including crystal field 
effects and application to real materials. One of the 
main bottlenecks in methods based on spectral moment 
expansions is the evaluation of high-order correlation 
functions. We find that including such correlations that 
are beyond two-particle type through approximate meth¬ 
ods such as CPA or lower order decomposition, can lead 
to spurious features at the chemical potential. We find 
the best benchmarks simply by neglecting correlations 
beyond two-particle. We conjecture that evaluation of 
the higher-order correlations through exact methods 
such as ligand held theory might be able to circumvent 
the issues mentioned aboveSSlS21_ ^re presently 

implementing such a procedure. This procedure will 
also enable us to treat the Hund’s coupling term in the 
rotationally invariant form rather than the simpler and 
approximate density-density type treated in the present 
work. We are also planning to extend our method to 
incorporate the off-diagonal hybridization, which at 
present is very difficult to handle by numerical exact 
methods like CTQMG. Apart from the benchmarks 
for model Hamiltonians in various parameter regimes, 
we have also carried out a GGA-I-DMFT(MO-IPT) 
study of the perovskite SrVOa, and compared our 
predicted photoemission spectrum with experiments 
and results from other methods. The agreement with 
experiments was found to be excellent. A full scale 
implementation of the method outlined here, with de¬ 
tailed instructions for installation and use may be found 
at http://www.institute.loni.org/lasigma/package/mo- 
ipt/ 


V. Acknowledgments 

We thank CSIR and DST (India) for research funding. 
Our simulations used an open source implementatiorpS 














17 


of the hybridization expansion continuous-time quan¬ 
tum Monte Carlo algorithnP^ and the ALP^^Hl libraries. 
This work is supported by NSF DMR-1237565 and NSF 
EPSCoR Cooperative Agreement EPS-1003897. Super¬ 
computer support is provided by the Louisiana Opti¬ 
cal Network Initiative (LONI) and HPC@LSU. NSV ac¬ 
knowledges an international travel fellowship award from 
lUSSTF. Dasari acknowledges the hospitality of the de¬ 
partment of Physics & Astronomy and the Center for 
Computation & Technology, at Louisiana State Univer¬ 
sity. 


VI. Appendix 


In this appendix, we provide the derivations of the un¬ 
known parameters appearing in the MO-IPT ansatz for 
the self energy (Eq. ©)• 

Derivation for Aq,: 


The spectral representation of the Q;*''-orbital Green’s 
function is given by 


G““(z) = f 

J —( 


D°‘°‘{e)de 


z — e 


(36) 


This can be Taylor expanded to obtain the Green func¬ 
tion in terms of spectral moments. 


J —oo ^ 


e e 


+...) 


^n+l 


n—0 


(37) 

where /j,„’s are the spectral moments. We can also repre¬ 
sent the Green function in terms of a continued fraction 
expansion and this is given by 


G““(z) = 


ai 


(38) 


1 + 


z + a^H- 


By comparing Eq. (38) with Eq. (37), we obtain the con¬ 


tinued fraction expansion coefficients in terms of spec¬ 
tral moments. Now we can calculate the spectral mo¬ 
ments exa ctly up to any order by using the following 

expressions^Sl 


a2 =-[(ea - fj-) + ^ Ua/3(n/3}j, (40) 

/3/(a) 


as = — 


- (drf 


CXCt , .CkCK 


dT MS 


(41) 


^ 2 “ 


— (Eq JV ^ ' ^ka 4 " 2(60, fl) ^ ( Uap{nji) 


+ EE , (42) 

/3^a 7^a 


as = 


(Eq, — fj,) Udjsinp) 


(43) 

For sufficiently large values of z, one can truncate the 
continued fraction expansion of the Green’s function 
(Eq. (38)) at the appropriate level and take the limit 
z —^ 00 . Up to the second order moment 


G““(z) = 


ai 


Q!2 


(44) 


After substituting the continued fraction expansion coef¬ 
ficients in Eq. (44), we find the self energy contribution 
to the Green’s function in the high frequency limit as 


E„(cc) 44^ C7„;3(n/3) 

/3/a 

X)^,7/q; UafiUaj {{nprij) — {'np){n^)) 
LO 


Ea(w) = > Uapirifs) H-- - - 

P=jta 

^ UapUa-y {{npn^) — {np){nj)) 

LO 

(45) 

In the high frequency limit the self energy ansatz reduces 
to the following form: 




OlOL 

n 




(n-p)-fold 


p-fold 



n = 0,1,2,.... ; 0< p < n 

The relation between the first few spectral moments and 
the continued fraction expansion coefficients is given by, 

= Mr = ({/«./!}) = 1, (39) 


E„(o;)= ^ U^0{np)+A^ ^ sg]. (46) 

/3/(a) /3/(a) 

It is easy to show that in the limit of high frequencies, 
^a ^/3 following fornP^, 

sS = —W(l-(^o/3))- (47) 

UJ 

Here ng^ is the Hartree-corrected impurity charge be¬ 
cause the propagators used in the second order pair bub¬ 
ble diagram are Hartree-corrected propagators. We ob- 
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tain the expression for Aa by substituting Eq. (471 in 
Eq. ([ 4 ^ and comparing with Eq. (451 as 


- (Mrf 

- (mD" ■ 


'^l3=jL{a) ^“7 ((^/S^t) ~ {'"^ 0 ) 

E/35^(a) t^a/sK/?) (1 - (no/ 3 )) 

Note that a two-particle correlation function is needed to 
find tIq,. 

Derivation for Bq,: 

The relation between the impurity Green’s function 
and the self energy in the atomic limit is, 


G„(a;) = 


1 


w~*" 01 — Cq — 


(48) 


where the self-energy, may be represented as a 

continued fraction: 


S„(a;) = + 0 — Ca — 


1 


(49) 


2 +^ 


2 + a4 • 


As a simple case we consider only two poles in the self 
energy. In principle we can keep all the poles of the self 
energy but the difficulty is that a pole of order n in¬ 
volves the (n -I- 1)*^ order correlation function. These 
functions are very hard to calculate without making ap¬ 
proximations. With the two pole approximation for the 
self energy, Eq. (|4^ reduces to the following fornP^ 


^a(^) — ^ -\- 


020:3 


a;+ -I- 03 -I- 04 


where 02 = —“ , 


03 = —- 


0T - ipr? 




(50) 


(51) 


(52) 


and 04 = — 


(53) 


In the atomic limit (V—>■ 0), the second order pair bubble 
diagram E^^j(a;) reduces to the following fornPUlll^ 


^(2), . G2^[(»^0/3)(1- (no/ 3 ))] 

= - 7]+TT, - 

UJ^ + 00 


(54) 


Here 00 is the pseudo-chemical potential. As mentioned 
earlier, we find this quantity by satisfying the Luttinger’s 
theorem or equivalently the Friedel’s sum rule. Now, the 
self energy ansatz becomes 

Eq, = 'y ] Uap{np) -\- 
/3/(a) 

^aE/3/(Q) Ulp[{nop)(l - (no/3))] 


UJ+ + 00 - Bq E/ 3 /(q) Gq^[(«0 /3)(1 - (no/ 3 ))] 


155) 


By comparing Eq. (551 with Eq. (49) we find the expres¬ 


sion for Bq in terms of spectral moments as, 
00 — (03 -|- 04 ) 


B„ = 


E,S/(a) ^a,s[(^0^)(l (nO/ 3 ))] 


(56) 


After substituting the spectral moments in Eq. (56 1 Be 
becomes. 


Bq = 


Mo T Cq 0 E/3^q GQ/3(n/3) 


E/35^q 'yZr):^a UapUa^Uar) [(n/3) (n-yU^) (n/3n-),n^)] 

(57) 




where 


= Ulp{nop){l - {nop)). (58) 

p^a 
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